library(quantreg)
library(affy)
data1<-read.table("frma_normalized1.txt",sep="\t",head=T)
data2<-read.table("frma_normalized2.txt",sep="\t",head=T)
data3<-read.table("frma_normalized3.txt",sep="\t",head=T)
data4<-read.table("frma_normalized4.txt",sep="\t",head=T)
data5<-read.table("frma_normalized5.txt",sep="\t",head=T)
data6<-read.table("frma_normalized6.txt",sep="\t",head=T)

Cen_1F<-apply(data1,1,mean)
Cen_2F<-apply(data2,1,mean)
Cen_3F<-apply(data3,1,mean)
Cen_4F<-apply(data4,1,mean)
Cen_5F<-apply(data5,1,mean)
Cen_6F<-apply(data6,1,mean)
CenF<-(Cen_1F+Cen_2F+Cen_3F+Cen_4F+Cen_5F+Cen_6F)/6
betweenF <- 20*((Cen_1F-CenF)^2+(Cen_2F-CenF)^2+(Cen_3F-CenF)^2+(Cen_4F-CenF)^2+(Cen_5F-CenF)^2+(Cen_6F-CenF)^2)/5
dataF<-as.matrix(cbind(data1[,1:20],data2[,1:20],data3[,1:20],data4[,1:20],data5[,1:20],data6[,1:20]))

acrossF<-CenF
for(i in 1:nrow(dataF)){
acrossF[i]<-0
}
for(i in 1:ncol(dataF)){
	if(i<21)
	{
		temp<-((dataF[,i]-Cen_1F)^2)/114
		acrossF<-acrossF + temp
	}
	if(i<41 &&i>20)
	{
		temp<-((dataF[,i]-Cen_2F)^2)/114
		acrossF<-acrossF + temp

	}
	if(i<61 && i>40)
	{
		temp<-((dataF[,i]-Cen_3F)^2)/114
		acrossF<-acrossF + temp

	}
	if(i<81 && i>60)
	{
		temp<-((dataF[,i]-Cen_4F)^2)/114
		acrossF<-acrossF + temp

	}
	if(i<101 && i>80)
	{
		temp<-((dataF[,i]-Cen_5F)^2)/114
		acrossF<-acrossF + temp

	}
	if(i<121 && i>100)
	{
		temp<-((dataF[,i]-Cen_6F)^2)/114
		acrossF<-acrossF + temp
	
	}
}
F<-betweenF/acrossF
data<-read.table("rma_normalizaed.txt",sep="\t",head=T)

Cen_1<-apply(data[,1:20],1,mean)
Cen_2<-apply(data[,21:40],1,mean)
Cen_3<-apply(data[,41:60],1,mean)
Cen_4<-apply(data[,61:80],1,mean)
Cen_5<-apply(data[,81:100],1,mean)
Cen_6<-apply(data[,101:120],1,mean)
Cen<-(Cen_1+Cen_2+Cen_3+Cen_4+Cen_5+Cen_6)/6
between <- 20*((Cen_1-Cen)^2+(Cen_2-Cen)^2+(Cen_3-Cen)^2+(Cen_4-Cen)^2+(Cen_5-Cen)^2+(Cen_6-Cen)^2)/5
across<-Cen
for(i in 1:nrow(data)){
across[i]<-0
}
for(i in 1:ncol(data)){
	if(i<21)
	{
		temp<-((data[,i]-Cen_1)^2)/114
		across<-across + temp
	}
	if(i<41 &&i>20)
	{
		temp<-((data[,i]-Cen_2)^2)/114
		across<-across + temp

	}
	if(i<61 && i>40)
	{
		temp<-((data[,i]-Cen_3)^2)/114
		across<-across + temp

	}
	if(i<81 && i>60)
	{
		temp<-((data[,i]-Cen_4)^2)/114
		across<-across + temp

	}
	if(i<101 && i>80)
	{
		temp<-((data[,i]-Cen_5)^2)/114
		across<-across + temp

	}
	if(i<121 && i>100)
	{
		temp<-((data[,i]-Cen_6)^2)/114
		across<-across + temp
	
	}
}
F1<-between/across


library(gplots)
#h1 <- hist(F1,breaks=seq(0,10,0.1),plot=FALSE)
#h2 <- hist(F,breaks=seq(0,10,0.1),plot=FALSE)
h1 <- hist(F1,breaks=seq(0,21,0.1),plot=FALSE)
h2 <- hist(F,breaks=seq(0,21,0.1),plot=FALSE)
hh <- cbind(h1$density,h2$density)
colnames(hh) <- c("RMA","fRMA")
mp <-barplot2(t(hh),beside=TRUE,col=c("red","blue"),legend=c("RMA","fRMA"),main="RMA vs fRMA F-statistics",xlim=c(0,400),ylab="Counts",plot.grid=TRUE)
arrows(x0=2.29,y0=0,y1=h2$density[-1],col="red",length=0.1)
box()

Top THREE QUARTILES
vec<-which(Cen>4.339362)
vecF<-which(Cen>4.577038)
comb_vec<-intersect(vec,vecF)
F<-F[comb_vec]
F1<-F1[comb_vec]
h1 <- hist(F1,breaks=seq(0,22516,0.1),plot=FALSE)
h <- hist(F,breaks=seq(0,22516,0.1),plot=FALSE)
h2 <- hist(F,breaks=seq(0,22516,0.1),plot=FALSE)
hh <- cbind(h1$density,h2$density)
colnames(hh) <- c("RMA","fRMA")
mp <-barplot2(t(hh),beside=TRUE,col=c("red","blue"),legend=c("RMA","fRMA"),main="RMA vs fRMA F-statistics",xlim=c(0,400),ylab="Counts",plot.grid=TRUE)




